Figure 8.31 shows the ACFs for 36 random numbers, 360 random numbers and 1,000 random numbers.
- Explain the differences among these figures. Do they all indicate that the data are white noise?
The figures show different critical values (blue dashed lines).
All figures indicate that the data are white noise.
- Why are the critical values at different distances from the mean of zero? Why are the autocorrelations different in each figure when they each refer to white noise?
A classic example of a non-stationary series is the daily closing IBM stock price series (data set
ibmclose). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.
The time plot shows the series “wandering around”, which is a typical indication of non-stationarity. Differencing the series should remove this feature.
ACF does not drop quickly to zero, moreover the value \(r_1\) is large and positive (almost 1 in this case). All these are signs of a non-stationary time series. Therefore it should be differenced to obtain a stationary series.
PACF value \(r_1\) is almost 1. All other values \(r_i, i>1\) are smaller than the critical value. These are signs of a non-stationary process that should be differenced in order to obtain a stationary series.
For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.
usnetelec
There is no need for a Box-Cox transformation in this case.
enplanements
visitors
For the
enplanementsdata, write down the differences you chose above using backshift operator notation.
A seasonal (lag 12) difference, followed by a first difference: \[ (1-B) (1-B^{12})\]
myts <- readxl::read_excel("data/retail.xlsx", skip=1)[,"A3349873A"] %>%
ts(frequency=12, start=c(1982,4))First we will try logs as before, followed by a seasonal difference.
That still looks a little non-stationary, so I will take another difference at lag 1.
That looks stationary now. So we used a log, a seasonal difference, and a first difference. Other time series from retail data set may require a different set of transformation/differencing choices.
Use R to simulate and plot some data from simple ARIMA models. a. Use the following R code to generate data from an AR(1) model with \(\phi_{1} = 0.6\) and \(\sigma^2=1\). The process starts with \(y_1=0\).
ar1 <- function(phi, n=100)
{
y <- ts(numeric(n))
e <- rnorm(n)
for(i in 2:n)
y[i] <- phi*y[i-1] + e[i]
return(y)
}
- Produce a time plot for the series. How does the plot change as you change \(\phi_1\)?
Some examples of changing \(\phi_1\)
- Write your own code to generate data from an MA(1) model with \(\theta_{1} = 0.6\) and \(\sigma^2=1\).
ma1 <- function(theta, n=100)
{
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y[i] <- theta*e[i-1] + e[i]
return(y)
}
- Produce a time plot for the series. How does the plot change as you change \(\theta_1\)?
- Generate data from an ARMA(1,1) model with \(\phi_{1} = 0.6\), \(\theta_{1} = 0.6\) and \(\sigma^2=1\).
arma11 <- function(phi, theta, n=100)
{
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100)
y[i] <- phi*y[i-1] + theta*e[i-1] + e[i]
return(y)
}
autoplot(arma11(0.6,0.6))
- Generate data from an AR(2) model with \(\phi_{1} =-0.8\), \(\phi_{2} = 0.3\) and \(\sigma^2=1\). (Note that these parameters will give a non-stationary series.)
ar2 <- function(phi1, phi2, n=100)
{
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 3:100)
y[i] <- phi1*y[i-1] + phi2*y[i-2] + e[i]
return(y)
}
autoplot(ar2(-0.8,0.3))
- Graph the latter two series and compare them.
See graphs above. The non-stationarity of the AR(2) process has led to increasing oscillations
Consider
wmurders, the number of women murdered each year (per 100,000 standard population) in the United States.
- By studying appropriate graphs of the series in R, find an appropriate ARIMA(\(p,d,q\)) model for these data.
The graphs suggest differencing of the data before applying ARMA models.
The ACF and PACF graphs suggest that the resulting time series might be an AR(2) or MA(2) due to the significant spike at lag 2 in both graphs. So either an ARIMA(0,1,2) or an ARIMA(2,1,0) appear appropriate. As we need to make a choice, I’ll use an ARIMA(0,1,2) in what follows.
- Should you include a constant in the model? Explain.
A constant would imply a drift in the original data which does not look correct, so we omit the constant.
- Write this model in terms of the backshift operator.
ARIMA(0,1,2): \((1-B)y_t = (1 + \theta_1B + \theta_2B^2) \varepsilon_{t}\).
- Fit the model using R and examine the residuals. Is the model satisfactory?
## Series: wmurders
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.066 0.371
## s.e. 0.126 0.164
##
## sigma^2 estimated as 0.0422: log likelihood=9.71
## AIC=-13.43 AICc=-12.95 BIC=-7.46
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)
## Q* = 9.77, df = 8, p-value = 0.28
##
## Model df: 2. Total lags used: 10
The model is satisfactory.
- Forecast three times ahead. Check your forecasts by hand to make sure that you know how they have been calculated.
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 2005 2.4585 2.1952 2.7217 2.0558 2.8611
## 2006 2.4771 2.1169 2.8373 1.9262 3.0280
## 2007 2.4771 1.9793 2.9749 1.7157 3.2385
To check the forecasts by hand we use the model written in terms of the backshift operator: \[ (1-B)y_t = (1 + \theta_1B + \theta_2B^2) \varepsilon_{t}, \] where \(y_t\) is the observed time series and the coefficients are \(\theta_1 = -0.066\) and \(\theta_2 = 0.371\). It can be rewritten as: \[ y_t = y_{t-1} + \varepsilon_t -0.066\varepsilon_{t-1} + 0.371\varepsilon_{t-2}. \] The last observation is 2.5894 and the last few residuals are -0.39266, -0.34376, 0.05024. So the next forecast is: \[ \hat{y}_{T+1|T} = 2.5894 + 0 -0.066 (0.0502) + 0.371 (-0.3438) = 2.4585. \] Similarly, \[ \hat{y}_{T+2|T} = 2.4585 + 0 + 0 + 0.371 (0.0502) = 2.4771 \] and \[ \hat{y}_{T+3|T} = 2.4771 + 0 + 0 + 0 = 2.4771. \]
- Create a plot of the series with forecasts and prediction intervals for the next three periods shown.
- Does
auto.arimagive the same model you have chosen? If not, which model do you think is better?
auto.arima picked a more complex model which assumes that the trend in the data will continue in the future. The extra level of differencing is hard to justify from the graphs, and it is unlikely that murder rates will fall indefinitely. So I’d probably prefer only one difference. We can force auto.arima to do that:
## Series: wmurders
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2
## -0.066 0.371
## s.e. 0.126 0.164
##
## sigma^2 estimated as 0.0422: log likelihood=9.71
## AIC=-13.43 AICc=-12.95 BIC=-7.46
Now auto.arima picks the same model that I did.
Consider
austa, the total international visitors to Australia (in millions) for the period 1980-2015.
- Use
auto.arima()to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.
## Series: austa
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## 0.301 0.173
## s.e. 0.165 0.039
##
## sigma^2 estimated as 0.0338: log likelihood=10.62
## AIC=-15.24 AICc=-14.46 BIC=-10.57
An ARIMA(0,1,1) with drift model is selected .
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1) with drift
## Q* = 2.3, df = 5, p-value = 0.81
##
## Model df: 2. Total lags used: 7
The residuals look like white noise.
- Plot forecasts from an ARIMA(0,1,1) model with no drift and compare these to part a. Remove the MA term and plot again.
Removing the drift term makes a big difference as the forecasts no longer have a linear trend. Removing the MA(1) term makes almost no difference.
- Plot forecasts from an ARIMA(2,1,3) model with drift. Remove the constant and see what happens.
Error in stats::arima(x = x, order = order, seasonal = seasonal, include.mean = include.mean, :
non-stationary AR part from CSS
The model fitting causes an error because the drift is really needed. What happens is that the AR terms try to handle the drift but end up going outside the stationarity region.
- Plot forecasts from an ARIMA(0,0,1) model with a constant. Remove the MA term and plot again.
A stationary model with a constant has long term forecasts equal to the mean (long term here meaning after the first horizon). A stationary model without a constant has long term forecasts equal to zero.
- Plot forecasts from an ARIMA(0,2,1) model with no constant.
The second level of differencing induces a linear trend in the forecasts, but the prediction intervals become quite wide.
For the
usgdpseries: a. if necessary, find a suitable Box-Cox transformation for the data;
I don’t think a Box-Cox transformation is required.
- fit a suitable ARIMA model to the transformed data using
auto.arima();
## Series: usgdp
## ARIMA(2,2,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## -0.123 0.311 -0.584 -0.367
## s.e. 0.287 0.087 0.300 0.286
##
## sigma^2 estimated as 1604: log likelihood=-1199.6
## AIC=2409.1 AICc=2409.4 BIC=2426.4
- try some other plausible models by experimenting with the orders chosen;
The second order differencing will induce a trend in the forecasts, which is required here, so I will look at changing only \(p\) and \(q\).
fit020 <- Arima(usgdp, order=c(0,2,0))
fit021 <- Arima(usgdp, order=c(0,2,1))
fit022 <- Arima(usgdp, order=c(0,2,2))
fit023 <- Arima(usgdp, order=c(0,2,3))
fit120 <- Arima(usgdp, order=c(1,2,0))
fit121 <- Arima(usgdp, order=c(1,2,1))
fit122 <- Arima(usgdp, order=c(1,2,2))
fit123 <- Arima(usgdp, order=c(1,2,3))
fit220 <- Arima(usgdp, order=c(2,2,0))
fit221 <- Arima(usgdp, order=c(2,2,1))
fit222 <- Arima(usgdp, order=c(2,2,2))
fit223 <- Arima(usgdp, order=c(2,2,3))
fit320 <- Arima(usgdp, order=c(3,2,0))
fit321 <- Arima(usgdp, order=c(3,2,1))
fit322 <- Arima(usgdp, order=c(3,2,2))
fit323 <- Arima(usgdp, order=c(3,2,3))
- choose what you think is the best model and check the residual diagnostics;
The best according to the AICc values is the ARIMA(2,2,1) model.
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,2,1)
## Q* = 9.36, df = 5, p-value = 0.095
##
## Model df: 3. Total lags used: 8
The residuals pass the Ljung-Box test, but the histogram looks like it has heavier tails than Gaussian.
- produce forecasts of your fitted model. Do the forecasts look reasonable?
These look reasonable.
- compare the results with what you would obtain using
ets()(with no transformation).
The ETS point forecasts are higher than the ARIMA forecasts, and the ETS forecast intervals are wider.
Consider
austourists, the quarterly number of international tourists to Australia for the period 1999–2010.
- Describe the time plot.
BoxCox.lambda function suggests using Box-Cox transformation with parameter \(\lambda=0.064\). This value is very close to \(0\) and moreover Box-Cox transformation with parameter \(0\) has another meaning: it is log transformation which ensures that the forecasted values transformed back to the original scale (using exponent transformation) will be positive. Since number of tourists visiting Australia cannot be negative it makes sense. Therefore we will just use logs of the original data.
The ACF and PACF graphs suggest seasonal differencing of the data before applying ARMA models.
The ACF and PACF graphs suggest that the resulting time series can be ARMA(0,1)(0,1) or ARMA(0,1)(1,0) or ARMA(1,0)(0,1) or ARMA(1,0)(1,0) models. Therefore the suggested ARIMA model for log(austourists) data are ARIMA(0,0,1)(0,1,1) or ARIMA(0,0,1)(1,1,0) or ARIMA(1,0,0)(0,1,1) or ARIMA(1,0,0)(1,1,0). A constant should be included in the model since the differenced time series is raised significantly above zero. The constant will lead to a linear trend in the forecasts. Such behaviour is strongly supported by the data.
fit1 = Arima(austourists, c(0,0,1), c(0,1,1), include.drift = TRUE, lambda = 0)
fit2 = Arima(austourists, c(0,0,1), c(1,1,0), include.drift = TRUE, lambda = 0)
fit3 = Arima(austourists, c(1,0,0), c(0,1,1), include.drift = TRUE, lambda = 0)
fit4 = Arima(austourists, c(1,0,0), c(1,1,0), include.drift = TRUE, lambda = 0)The best model out of these four can be chosen using AICc:
## [1] -167.84 -156.36 -170.40 -159.63
The lowest AICc is for model ARIMA(1,0,0)(0,1,1) with drift, therefore it is the preferred model for forecasting.
## Series: austourists
## ARIMA(1,0,0)(2,1,0)[4] with drift
## Box Cox transformation: lambda= 0
##
## Coefficients:
## ar1 sar1 sar2 drift
## 0.433 -0.710 -0.237 0.013
## s.e. 0.118 0.126 0.134 0.002
##
## sigma^2 estimated as 0.00415: log likelihood=85.65
## AIC=-161.3 AICc=-160.27 BIC=-150.51
auto.arima finds the same model according to AICc.
The best model written using backshift operator and without it
Using backshift operator the model can be written as: \((1 - \phi_1 B) (1 - B^4) y_t = c + (1 + \theta_1 B^4) \epsilon_t\)
Without the backshift operator the model can be written as: \(y_t - \phi_1 y_{t-1} - y_{t-4} + \phi_1 y_{t-5} = c + \epsilon_t + \theta_1 \epsilon_{t-4}\) or: \(y_t = \phi_1 y_{t-1} + y_{t-4} - \phi_1 y_{t-5} + c + \epsilon_t + \theta_1 \epsilon_{t-4}\)
The last form allows forecasting value \(y_t\) assuming that values \(y_{t-1}\), \(y_{t-4}\), \(y_{t-5}\) and \(\epsilon_{t-4}\) are known. Value \(\epsilon_t\), although, should be estimated as zero.
- What can you learn from the ACF graph?
- What can you learn from the PACF graph?
- Produce plots of the seasonally differenced data \((1 - B^{4})Y_{t}\). What model do these graphs suggest?
- Does
auto.arima()give the same model that you chose? If not, which model do you think is better?
- Write the model in terms of the backshift operator, then without using the backshift operator.
Consider
usmelec, the total net generation of electricity (in billion kilowatt hours) by the U.S. electric industry (monthly for the period January 1973 – June 2013). In general there are two peaks per year: in mid-summer and mid-winter.
- Examine the 12-month moving average of this series to see what kind of trend is involved.
- Do the data need transforming? If so, find a suitable transformation.
- Are the data stationary? If not, find an appropriate differencing which yields stationary data.
- Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?
- Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
- Forecast the next 15 years of electricity generation by the U.S. electric industry. Get the latest figures from the EIA to check the accuracy of your forecasts.
- Eventually, the prediction intervals are so wide that the forecasts are not particularly useful. How many years of forecasts do you think are sufficiently accurate to be usable?
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1973 160.22 143.54 148.16 139.59 147.40 161.24 173.73 177.37 156.88 154.20
## 1974 157.56 142.75 150.34 142.31 153.81 156.44 178.25 174.12 152.47 152.20
## 1975 164.62 147.35 155.76 146.50 153.53 162.72 177.06 179.93 155.44 155.19
## 1976 178.61 156.97 164.47 153.47 157.66 173.67 186.69 186.64 165.24 164.01
## 1977 196.66 162.95 169.44 157.12 169.60 181.03 199.17 196.36 176.50 166.65
## 1978 198.11 173.75 173.46 160.01 175.55 188.59 202.95 206.66 185.80 176.01
## 1979 209.99 186.59 183.15 170.26 178.41 186.98 202.52 205.10 180.97 179.95
## 1980 200.30 188.96 187.75 169.02 176.07 189.75 217.06 215.63 191.70 178.76
## 1981 206.76 179.86 185.83 172.84 178.14 203.02 220.66 210.64 187.05 181.56
## 1982 209.69 180.55 187.97 172.88 177.48 186.45 210.87 205.89 180.88 173.17
## 1983 195.87 172.72 182.77 170.67 174.72 191.37 220.45 230.19 195.82 183.14
## 1984 216.92 189.81 200.39 181.38 192.55 209.97 221.53 229.53 195.41 191.14
## 1985 228.15 198.49 195.25 185.17 197.12 205.68 227.00 226.29 202.71 195.00
## 1986 217.76 192.58 197.12 186.37 197.65 215.33 242.95 225.40 206.91 197.96
## 1987 223.04 194.28 202.13 189.79 206.41 225.91 248.20 247.88 213.22 203.22
## 1988 238.19 217.18 214.29 196.30 208.70 233.07 257.74 267.93 220.39 210.81
## 1989 246.77 233.78 241.95 222.87 234.82 250.92 273.46 275.42 242.79 235.67
## 1990 255.19 229.50 244.76 229.77 241.77 268.99 287.45 290.65 258.36 244.37
## 1991 269.21 228.29 240.56 227.80 254.87 269.10 294.65 292.16 256.03 245.38
## 1992 267.77 239.51 247.73 233.41 242.41 261.08 293.62 281.93 259.93 244.99
## 1993 271.02 248.01 261.25 234.69 244.33 275.36 312.23 311.45 264.03 250.55
## 1994 289.77 249.17 258.00 240.64 252.75 294.16 311.26 307.61 266.26 256.53
## 1995 279.77 252.31 261.34 244.74 264.29 286.26 330.42 345.78 277.57 263.98
## 1996 296.92 270.69 275.02 251.61 282.27 302.72 327.71 329.29 283.15 271.45
## 1997 300.57 258.13 272.26 258.28 272.91 299.09 344.52 332.90 301.06 284.97
## 1998 295.26 260.59 286.88 261.23 299.64 328.90 361.94 357.37 318.92 284.45
## 1999 315.81 274.82 298.14 280.72 300.10 328.92 376.54 364.03 305.52 283.94
## 2000 327.99 294.17 301.58 285.58 322.95 339.05 356.53 368.67 312.45 289.45
## 2001 332.49 282.94 300.71 278.08 300.49 327.69 357.61 370.53 306.93 294.73
## 2002 319.94 281.83 302.55 289.85 307.68 341.02 381.54 374.59 331.28 307.06
## 2003 341.99 299.25 304.32 285.76 307.55 328.69 374.40 381.82 323.14 306.74
## 2004 346.55 314.28 308.81 290.56 327.38 345.08 377.33 368.44 335.62 312.45
## 2005 343.12 298.50 317.46 289.56 315.06 363.67 402.27 404.94 350.22 316.40
## 2006 328.66 307.33 318.73 297.86 330.62 364.26 410.42 407.76 332.06 321.57
## 2007 353.53 323.23 320.47 303.13 330.20 362.75 393.23 421.80 355.39 332.62
## 2008 363.00 325.11 324.63 305.87 325.25 373.11 402.90 388.99 338.06 318.55
## 2009 354.99 300.89 310.60 289.54 311.31 347.66 372.54 381.22 327.40 307.04
## 2010 360.96 319.74 312.17 287.80 327.94 375.76 409.73 408.88 346.05 307.92
## 2011 363.11 313.29 318.71 302.40 323.63 367.73 418.69 406.54 337.96 308.73
## 2012 340.92 310.15 309.04 295.94 337.53 361.51 416.51 396.11 334.74 312.16
## 2013 348.64 309.60 325.37 298.26 322.12 356.40
## Nov Dec
## 1973 148.14 153.60
## 1974 150.07 160.01
## 1975 153.03 169.63
## 1976 169.35 184.14
## 1977 167.39 184.59
## 1978 176.39 192.10
## 1979 177.77 188.97
## 1980 178.77 195.85
## 1981 175.79 195.83
## 1982 173.60 184.96
## 1983 183.17 212.56
## 1984 190.60 200.23
## 1985 192.65 219.49
## 1986 196.65 213.79
## 1987 200.48 220.74
## 1988 209.81 232.99
## 1989 234.25 274.45
## 1990 231.24 255.77
## 1991 241.66 254.09
## 1992 244.16 267.34
## 1993 252.10 272.17
## 1994 251.87 269.52
## 1995 261.48 285.55
## 1996 269.12 284.25
## 1997 271.12 296.36
## 1998 267.14 297.99
## 1999 271.10 295.17
## 2000 284.67 319.01
## 2001 278.93 305.50
## 2002 296.29 324.83
## 2003 297.87 331.68
## 2004 302.10 341.95
## 2005 306.12 348.10
## 2006 309.16 336.28
## 2007 314.10 346.29
## 2008 310.05 343.90
## 2009 296.63 350.51
## 2010 306.01 362.12
## 2011 304.12 335.75
## 2012 305.55 334.33
## 2013
The data show strong up-trend which might be considered as linear in long run.
The seasonal pattern which size increases together with the level of the trend strongly suggests using Box-Cox transformation before forecasting.
BoxCox.lambda function suggests using Box-Cox transformation with parameter \(\lambda=-0.57\).
The data is not stationary (it contains trend), therefore differencing is required to make it stationary for ARMA forecasting.
ACF and PACF graphs suggest that the data is stationary. Suggested models are: ARIMA(1,0,0)(0,1,1) ARIMA(2,0,0)(0,1,1) ARIMA(3,0,0)(0,1,1)
fit1 = Arima(usmelec, c(1,0,0), c(0,1,1), lambda = lambda)
fit2 = Arima(usmelec, c(2,0,0), c(0,1,1), lambda = lambda)
fit3 = Arima(usmelec, c(3,0,0), c(0,1,1), lambda = lambda)## [1] -5002.6 -5032.7 -5058.1
The lowest AICc is for model ARIMA(3,0,0)(0,1,1).
##
## Box-Ljung test
##
## data: fit3$residuals
## X-squared = 11.2, df = 1, p-value = 8e-04
ACF and PACF graphs as well as Ljung-Box test suggest that the residuals are not white noise. Looking at the residuals we can suggest spikes at lags 3 and 4 at ACF graph can be removed by increasing MA value in ARIMA model:
fit4 = Arima(usmelec, c(3,0,1), c(0,1,1), lambda = lambda)
fit5 = Arima(usmelec, c(3,0,2), c(0,1,1), lambda = lambda)
fit6 = Arima(usmelec, c(3,0,3), c(0,1,1), lambda = lambda)## [1] -5089.8 -5087.0 -5085.8
It appears that currently model ARIMA(3,0,1)(0,1,1) is with the lowest AICc.
The residuals for this model are:
Which are good except for the spikes at ACF and PACF graphs at lag 15. We will try to rectify it by increasing AR and MA values for seasonal and non-seasonal ARIMA parts:
fit7 = Arima(usmelec, c(3,0,1), c(1,1,1), lambda = lambda)
fit8 = Arima(usmelec, c(4,0,1), c(0,1,1), lambda = lambda)
fit9 = Arima(usmelec, c(4,0,1), c(1,1,1), lambda = lambda)## [1] -5089.1 -5087.8 -5087.2
All received AICc values are greater than for the best found model ARIMA(3,0,1)(0,1,1). Therefore ARIMA(3,0,1)(0,1,1) should be used for forecasting. The residuals although appear to be correlated for lag 15:
##
## Box-Ljung test
##
## data: fit4$residuals
## X-squared = 27.8, df = 11, p-value = 0.0034
It means that the forecasting intervals, calculated using this model can be incorrect.
# New data taken from http://www.eia.gov/totalenergy/data/monthly/#electricity website
newData = c(306009.629,362119.081,362871.893,313126.607,318709.941,302400.724,323628.238,
367727.015,418692.755,406511.315,337931.318,308698.504,304102.155,335740.463,
339528.259,309389.42,309090.569,295228.216,336517.631,360826.151,414639.671,
395699.566,334584.877,311650.774,305975.035,334635.089,348489.737,309435.119,
325301.486,298073.696,321833.617,356224.262,393799.21,383968.062,340293.27,
314682.729,313751.509,352356.534,377019.193,323662.011,331595.192,296766.045,
323730.672,357419.408,384838.868,383494.071,338975.899,313971.928,317176.051)/1000
newTs = ts(newData, start = c(2010, 11), frequency = 12)
fcast = forecast(fit4, h = 15*12)
plot(fcast)
lines(newTs, col = "red", lwd = 1.5)## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.13594 7.2298 5.3642 -0.014499 2.0398 0.59571 -0.030083
## Test set 1.05573 11.0216 8.5795 0.410227 2.4113 0.95278 0.593354
## Theil's U
## Training set NA
## Test set 0.37767
The forecasts are good. Probably quite a few years of forecasts can be correct. Although the forecasting intervals are expanding very quickly and therefore the forecasts become less and less useful when the forecasting horizon increases.
For the
mcopperdata:
- if necessary, find a suitable Box-Cox transformation for the data;
- fit a suitable ARIMA model to the transformed data using
auto.arima();
- try some other plausible models by experimenting with the orders chosen;
- choose what you think is the best model and check the residual diagnostics;
- produce forecasts of your fitted model. Do the forecasts look reasonable?
- compare the results with what you would obtain using
ets()(with no transformation).
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1960 255.2 259.7 249.3 258.0 244.3 246.8 250.6 241.3 231.0 218.7
## 1961 216.6 220.1 221.7 225.6 238.6 232.8 226.1 227.2 225.8 225.0
## 1962 226.8 231.3 231.1 230.6 230.5 230.4 230.4 230.4 230.4 230.6
## 1963 230.6 230.6 230.6 230.6 230.6 230.6 230.6 230.6 230.6 230.6
## 1964 234.2 247.8 266.1 307.7 295.8 288.8 305.4 356.9 415.0 485.1
## 1965 356.7 419.2 440.6 480.5 490.8 466.1 404.0 431.6 473.5 500.1
## 1966 599.0 668.7 668.7 679.9 592.8 604.8 559.5 426.4 402.4 454.9
## 1967 443.7 435.4 391.8 354.9 369.5 362.3 355.9 372.7 378.3 406.3
## 1968 586.4 715.8 708.0 522.4 456.5 473.0 439.9 439.9 462.3 449.7
## 1969 523.8 535.0 532.4 578.1 579.6 617.0 605.4 669.0 657.6 647.4
## 1970 677.2 690.1 730.4 725.2 665.9 607.0 567.8 527.5 519.3 475.8
## 1971 420.7 424.9 476.5 521.2 464.1 447.2 464.4 451.0 427.5 417.7
## 1972 418.9 426.9 442.2 433.5 423.3 412.4 423.2 427.0 434.2 428.7
## 1973 474.7 511.8 610.2 638.8 613.0 678.4 795.4 843.5 799.6 849.6
## 1974 913.2 1006.5 1172.0 1267.7 1190.6 1020.1 802.8 767.9 630.1 599.2
## 1975 512.4 528.7 554.5 560.5 539.7 522.5 559.3 603.7 580.2 573.1
## 1976 587.8 601.6 683.7 817.9 836.0 878.0 922.1 862.2 844.9 784.8
## 1977 814.8 833.4 881.0 830.4 797.9 763.1 724.4 665.8 685.6 682.9
## 1978 651.2 627.0 657.6 694.4 716.0 725.0 705.4 734.6 735.8 750.1
## 1979 827.1 969.7 1005.5 1012.3 935.3 888.8 802.4 883.2 953.6 967.6
## 1980 1147.9 1220.2 1045.5 937.1 888.9 858.5 916.7 877.9 857.6 846.2
## 1981 777.4 784.9 814.3 837.0 834.0 861.0 897.3 981.4 942.2 904.9
## 1982 853.6 863.8 836.3 858.5 843.6 740.4 829.9 841.0 832.7 861.4
## 1983 997.3 1074.3 1071.9 1090.0 1122.8 1098.6 1115.6 1091.0 1041.0 958.4
## 1984 978.0 991.5 1030.9 1078.1 1022.3 991.2 1007.9 1018.2 1029.3 1043.5
## 1985 1205.5 1269.8 1234.7 1212.7 1225.3 1118.0 1067.5 1025.7 1001.4 973.8
## 1986 995.4 982.8 984.4 957.3 932.3 937.0 891.7 876.5 916.1 922.9
## 1987 893.8 902.7 920.1 909.3 911.9 964.6 1052.8 1097.5 1100.7 1182.7
## 1988 1476.8 1324.3 1286.3 1216.5 1307.1 1428.8 1297.5 1296.2 1445.7 1689.4
## 1989 1912.6 1765.0 1904.1 1832.3 1679.1 1638.8 1539.1 1731.5 1834.9 1801.8
## 1990 1432.7 1390.9 1615.7 1639.6 1633.7 1510.5 1529.6 1554.4 1611.6 1409.6
## 1991 1265.1 1246.5 1322.5 1412.3 1336.8 1344.8 1353.8 1325.6 1346.1 1371.2
## 1992 1182.2 1240.5 1291.7 1260.9 1224.6 1239.1 1313.8 1297.2 1307.0 1360.3
## 1993 1472.3 1536.7 1472.3 1261.9 1159.0 1228.5 1287.5 1305.0 1219.7 1108.7
## 1994 1208.4 1261.0 1284.0 1268.1 1430.5 1549.9 1589.5 1560.0 1601.0 1586.2
## 1995 1910.5 1830.3 1826.8 1805.7 1745.9 1876.8 1927.4 1936.6 1870.3 1782.7
## 1996 1708.7 1650.2 1676.5 1713.5 1757.3 1407.7 1276.6 1294.9 1244.3 1236.2
## 1997 1469.2 1480.3 1506.4 1466.8 1540.0 1588.0 1466.7 1403.6 1315.4 1256.4
## 1998 1032.1 1014.4 1051.4 1078.1 1057.8 1005.6 1004.2 991.8 979.1 935.6
## 1999 866.7 866.5 849.6 910.2 935.6 891.3 1041.2 1025.5 1077.3 1040.1
## 2000 1124.2 1124.7 1100.6 1059.8 1183.3 1161.6 1192.3 1245.1 1365.3 1308.1
## 2001 1209.8 1214.9 1202.8 1159.3 1178.9 1147.5 1078.5 1019.0 974.4 948.5
## 2002 1049.9 1097.4 1127.9 1101.4 1092.9 1110.0 1022.2 962.4 950.0 952.5
## 2003 1018.8 1046.7 1047.4 1007.9 1015.5 1015.4 1054.0 1103.4 1109.1 1143.9
## 2004 1329.0 1477.5 1647.3 1635.6 1529.5 1469.7 1523.4 1563.1 1614.9 1666.7
## 2005 1687.8 1723.8 1773.0 1790.0 1751.7 1938.0 2063.8 2115.8 2133.3 2301.1
## 2006 2677.7 2851.6 2927.3 3610.9 4306.0 3897.7 4170.5 4059.8 4032.9 3998.6
## Nov Dec
## 1960 222.8 227.5
## 1961 225.7 226.3
## 1962 230.4 230.5
## 1963 230.6 232.2
## 1964 500.1 453.3
## 1965 523.8 541.4
## 1966 464.4 433.4
## 1967 515.3 551.3
## 1968 458.0 493.4
## 1969 679.4 707.8
## 1970 451.9 435.4
## 1971 406.4 410.7
## 1972 427.8 435.7
## 1973 950.4 959.9
## 1974 608.2 553.2
## 1975 575.1 568.9
## 1976 780.3 767.3
## 1977 650.2 679.1
## 1978 749.3 771.7
## 1979 978.3 1005.6
## 1980 839.5 800.4
## 1981 867.7 869.3
## 1982 884.3 911.4
## 1983 940.4 986.7
## 1984 1085.1 1113.5
## 1985 950.9 962.4
## 1986 915.0 927.6
## 1987 1419.9 1566.6
## 1988 1826.2 1915.3
## 1989 1647.3 1514.1
## 1990 1315.9 1292.7
## 1991 1336.8 1216.4
## 1992 1413.2 1422.5
## 1993 1100.3 1156.4
## 1994 1763.6 1914.4
## 1995 1904.8 1898.7
## 1996 1343.7 1366.9
## 1997 1135.0 1061.5
## 1998 946.9 881.9
## 1999 1065.3 1093.7
## 2000 1258.9 1264.4
## 2001 994.2 1020.8
## 2002 1006.4 1005.7
## 2003 1216.0 1258.5
## 2004 1678.5 1631.1
## 2005 2461.6 2621.7
## 2006 3675.9 3435.8
BoxCox.lambda function suggests using Box-Cox transformation with parameter \(\lambda=0.19\).
## Series: mcopper
## ARIMA(0,1,1)
## Box Cox transformation: lambda= 0.1919
##
## Coefficients:
## ma1
## 0.372
## s.e. 0.039
##
## sigma^2 estimated as 0.05: log likelihood=45.05
## AIC=-86.1 AICc=-86.08 BIC=-77.43
Other models which can be tried are ARIMA(2,1,0) and ARIMA(7,1,0):
## [1] -86.075 -82.512 -80.075
The lowest AICc is for model ARIMA(0,1,1) ( the model which was proposed by auto.arima).
##
## Box-Ljung test
##
## data: fit$residuals
## X-squared = 9.48, df = 11, p-value = 0.58
The residuals look normal at the histogram
No obvious patterns are visible on the graph of the residuals
The residuals are white noise
Forecasts look differently, but taking into account that the forecasting intervals are very wide, both forecasts provide similar results.
Choose one of the following seasonal time series:
hsales,auscafe,qauselec,qcement,qgas.
- Do the data need transforming? If so, find a suitable transformation.
- Are the data stationary? If not, find an appropriate differencing which yields stationary data.
- Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?
- Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
- Forecast the next 24 months of data using your preferred model.
- Compare the forecasts obtained using
ets().
For the same time series you used in the previous exercise, try using a non-seasonal model applied to the seasonally adjusted data obtained from STL. The
stlf()function will make the calculations easy (withmethod="arima"). Compare the forecasts with those obtained in the previous exercise. Which do you think is the best approach?
For your retail time series (Exercise 5 above):
- develop an appropriate seasonal ARIMA model;
- compare the forecasts with those you obtained in earlier chapters;
- Obtain up-to-date retail data from the ABS website (Cat 8501.0, Table 11), and compare your forecasts with the actual numbers. How good were the forecasts from the various models?
- Produce a time plot of the sheep population of England and Wales from 1867–1939 (data set
sheep).
- Assume you decide to fit the following model: \[ y_t = y_{t-1} + \phi_1(y_{t-1}-y_{t-2}) + \phi_2(y_{t-2}-y_{t-3}) + \phi_3(y_{t-3}-y_{t-4}) + e_t, \] where \(e_t\) is a white noise series. What sort of ARIMA model is this (i.e., what are \(p\), \(d\), and \(q\))?
ARIMA(3,1,0). \(p=3\), \(d=1\), \(q=0\).
- By examining the ACF and PACF of the differenced data, explain why this model is appropriate.
The last significant spike in the PACF is at lag 3, so \(p=3\) is a good choice. Hence an ARIMA(3,1,0) model.
- The last five values of the series are given below:
Year 1935 1936 1937 1938 1939 Millions of sheep 1648 1665 1627 1791 1797 The estimated parameters are \(\phi_1 = 0.42\), \(\phi_2 = -0.20\), and \(\phi_3 = -0.30\). Without using the
forecastfunction, calculate forecasts for the next three years (1940–1942).
From the equation above \[\begin{align*} \hat{y}_{T+1|T} &= y_{T} + \phi_1(y_{T}-y_{T-1}) + \phi_2(y_{T-1}-y_{T-2}) + \phi_3(y_{T-2}-y_{T-3})\\ &= 1797 + 0.42 (1797-1791) + -0.20 (1791-1627) + -0.30 (1627-1665) \\ &= 1778.00 \\ \hat{y}_{T+2|T} &= \hat{y}_{T+1|T} + \phi_1(\hat{y}_{T+1|T}-y_{T}) + \phi_2(y_{T}-y_{T-1}) + \phi_3(y_{T-1}-y_{T-2})\\ &= 1777.99577 + 0.42 (1777.99577 -1797) + -0.20 (1797-1791) + -0.30 (1791-1627) \\ &= 1718.87 \\ \hat{y}_{T+3|T} &= \hat{y}_{T+2|T} + \phi_1(\hat{y}_{T+2|T}-\hat{y}_{T+1|T}) + \phi_2(\hat{y}_{T+1|T}-y_{T}) + \phi_3(y_{T}-y_{T-1})\\ &= 1718.86855 + 0.42 (1718.86855 - 1777.99577) + -0.20 (1777.99577-1797) + -0.30 (1797-1791) \\ &= 1695.99 \end{align*}\]
- Now fit the model in R and obtain the forecasts using
forecast. How are they different from yours? Why?
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1940 1778.0 1687.5 1868.5 1639.5 1916.5
## 1941 1718.9 1561.5 1876.2 1478.3 1959.5
## 1942 1696.0 1494.2 1897.8 1387.3 2004.7
They are the same. If you had different values, it is most likely due to rounding errors.
- Plot the annual bituminous coal production in the United States from 1920 to 1968 (data set ).
- You decide to fit the following model to the series: \[y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \phi_3 y_{t-3} + \phi_4 y_{t-4} + e_t\] where \(y_t\) is the coal production in year \(t\) and \(e_t\) is a white noise series. What sort of ARIMA model is this (i.e., what are \(p\), \(d\), and \(q\))?
ARIMA(4,0,0). \(p=4\), \(d=0\), \(q=0\).
- Explain why this model was chosen using the ACF and PACF.
The PACF has a signifcant spike at lag 4, but none at higher lags.
- The last five values of the series are given below.
Year 1964 1965 1966 1967 1968 Millions of tons 467 512 534 552 545
The estimated parameters are \(c = 162.00\), \(\phi_1 = 0.83\), \(\phi_2 = -0.34\), \(\phi_3 = 0.55\), and \(\phi_4 = -0.38\). Without using the
forecastfunction, calculate forecasts for the next three years (1969–1971).
\[\begin{align*} \hat{y}_{T+1|T} &= c + \phi_1 y_{T} + \phi_2 y_{T-1} + \phi_3 y_{T-2} + \phi_4 y_{T-3} \\ &= 162.00 + 0.83 *545 -0.34 *552 + 0.55 *534 -0.38 *512 \\ &= 527.63 \\ \hat{y}_{T+2|T} &= c + \phi_1 \hat{y}_{T+1|T} + \phi_2 y_{T} + \phi_3 y_{T-1} + \phi_4 y_{T-2} \\ &= 162.00 + 0.83* 527.63 -0.34 *545 + 0.55 *552 -0.38* 534 \\ &= 517.19 \\ \hat{y}_{T+3|T} &= c + \phi_1 \hat{y}_{T+2|T} + \phi_2 \hat{y}_{T+1|T} + \phi_3 y_{T} + \phi_4 y_{T-1} \\ &= 162.00 + 0.83* 517.19 -0.34 * 527.63 + 0.55 *545 -0.38* 552 \\ &= 503.81 \\ \end{align*}\]
- Now fit the model in R and obtain the forecasts from the same model. How are they different from yours? Why?
## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 1969 527.63 459.88 595.38 424.02 631.24
## 1970 517.19 429.00 605.38 382.32 652.07
## 1971 503.81 412.48 595.13 364.13 643.48
Again, any differences between these values and those computed “by hand” are probably due to rounding errors.